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Abstract 

A parameter estimation method is devised for a slow-fast stochastic dynamical system, where 
often only the slow component is observable. By using the observations only on the slow compo- 
nent, the system parameters are estimated by working on the slow system on the random slow 
manifold. This offers a benefit of dimension reduction in quantifying parameters in stochastic 
dynamical systems. An example is presented to illustrate this method, and verify that the pa- 
rameter estimator based on the lower dimensional, reduced slow system is a good approximation 
of the parameter estimator for original slow-fast stochastic dynamical system. 

Mathematics Subject Classifications (2010) Primary 60H30; Secondary 60H10; 37D10. 

Keywords Parameter estimation; Slow-fast system; Random slow manifold; Quantifying uncer- 
tainty; Numerical optimization 



1 Introduction 

Invariant manifolds provide geometric structures for understanding dynamical behavior of nonlinear 
systems under uncertainty. Some systems evolve on fast and slow time scales, and may be modeled 
by coupled singularly perturbed stochastic ordinary differential equations (SDEs). A slow- fast 
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stochastic system may have a special invariant manifold called a random slow manifold that capture 
the slow dynamics. 

We consider a stochastic slow-fast system 



X = Ax + f{x, y) 



x(0) = xo G M' 



n 



(1) 



ij = -By + -g{x, y) + -^Wt, y(0) = yo G 




m 



(2) 



where A and B are matrices, e is a small positive parameter measuring slow and fast scale sep- 
aration, / and g are nonlinear Lipschitz continuous functions with Lipschitz constant Lf and Lg 
respectively, o" is a noise intensity constant, and {Wt : t E M} is a two-sided M™- valued Wiener 
process (i.e., Brownian motion) on a probability space (SI, J-", P). Under a gap condition and for e 
sufficiently small, there exists a random slow manifold y = h^{^,u}), uj £ 0,, as in [l?>, 5], for slow- 
fast stochastic system (l)-(2). When the nonlinearities f,g are only locally Lipschitz continuous 
but the system has a random absorbing set (e.g., in mean-square norm), we conduct a cut-off of 
the original system. The new system will have a random slow manifold which captures the original 
system's slow dynamics. 

The random slow manifold is the graph of a random nonlinear mapping h^{^, uj) = ar]^{uj)+ hf{^., oj), 
with /i^(^,a;) determined by a Lyapunov- Perron integral equation [5], 



with rj'^{9t(jj) = I-oo^^^^ dWs and r]^{Lo) = J-oo^ ^^dWs- The random slow manifold 

exponentially attracts other solution orbits. We will find an analytically approximated random 
slow manifold for sufficiently small e, in terms of an asymptotic expansion in e, as in [14, 15]. This 
slow manifold may also be numerically computed as in [7]. By restricting to the slow manifold, we 
obtain a lower dimensional reduced system of the original slow-fast system (l)-(2), for e sufficiently 
small 



where 9t and ips are defined in the next section. 

If the original slow- fast system (l)-(2) contains unknown system parameters, but only the slow 
component x is observable, we conduct parameter estimation using the slow system (3). Since 
the slow system is lower dimensional than the original system, this method offers an advantage in 
computational cost, in addition to the benefit of using only observations on slow variables. 

This paper is arranged as follows. In the next section, we obtain an approximated random slow 
manifold and thus the random slow system. Then in Section 3, we present a method for parameter 
estimation on the slow manifold. Finally, a simple example is presented in Section 4 to illustrate 
our method. 





n 



(3) 
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2 Random slow manifold and its approximation 

By a random transformation 

(j):=V.C.,.,-,)=(^_;^.(^)). (4) 

we convert the SDE system (l)-(2) to the foUowing system with random coefficients 

X{t) = AX{t) + f{X{t), Y{t) + ar,'{dtuj)), (5) 

Y{t) = ^BY{t) + ^g{Xit),Y{t) + ar]'i9tLo)), (6) 

where ri'^{uj) = J^^o dWg is the stationary solution of linear system dy'^ = ^y^dt + -^dWt- 
Here 0^ : 17 — > is the Wiener shift imphcitly defined by Ws{9tuj) = Wt+s{oj) — Wt{u}). Note that 

Define a mapping (between random samples) ip;; : Q, ^ Q implicitly by Wt{'ipe^) = -^Wte{^)- 
Then -^Wte{(^) is also a Wiener process with the same distribution as Wt{oj). Moreover, r}'^{9te(^) 

and rj^ijjj) are identically distributed with r]{9tip£u) = e^^''~^^ dWsii^e^) and ri{'ilj^uj) = e~^'^ dWs{'4'e'^)i 
respectively. 

By a time change r = t/e and using the fact that r]^{9T-e^^) and r}{9T-ij)f,Lo) are identically dis- 
tributed, the system (5)- (6) is reformulated as 

X' = e[AX + f{X,Y + ar^{9r^e^))l (7) 
Y' = BY + g{X,Y + a'n{9rij,uj)), (8) 

where ' = ^• 

We make the following two hypotheses. 

HI: There are positive constants a, /3 and K, such that for every x € M" and y € M"*, the 
following exponential estimates hold: 

|e^*x|M. < ire"*|x|M", t < 0; |e^*y|R™ < Ke-%\w.., t > 0. 
H2: (3 > KLg. 

Then there exists a random slow manifold M^(uj) = {{S^,h^{^,uj)) : ^ € M"} for the random 
system (5)-(6), with /i^ being expressed as follows [5, 12], 

h'{^,u) = - f e-T^g{X{s,uj,i),Y{s,,o,i) + arj'{9sUj))ds. 

^ J —CO 
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We can get a small e approximation for h'^. Start with the expansion Y = Yq + eYi + 0{e^) and the 
integral expression X = X{0) + £ f^[AX + f{X, Y + ar]{6ripe^))] dr in the system (7)-(8). Using the 
Taylor expansions of f{X, Y + ar]{9ripe^)) and g{X, Y + a7]{9rtpe^)) at point (xq, Yq + ar]{6T-ips^)) 
and denoting ^ = X{0), we obtain 

h%^,u;) = h^it^) + ehlituj) + O(e'), (9) 
for e sufficiently small, with 

m,uj)= r e~'''g{^,Yo{s)+ar,{9sA^))ds, (10) 



oo 



and 

rO 



m,u^) = / e-^^|(7..(e,lo(s)+^T7?(Me^))[Ase+ / f{^,Yo{r)+ar|{er^l^,u))dr] 

J-oo ^ Jo 

+gy{^,Yo{s) + a7]i9s^eUj))Yiis)}ds, (11) 
where Yq^I) and Yi{t) satisfy the following random differential equations, respectively 



Y^{t) = BYoir) + g{C, Yo{t) + ari{9rAuj)), 
ro(0) = /i^rf(0,a;), 



(12) 



and 



■y/(r) = [B + gyi^Yoir) + ar]{9^il;eCo))]Yi{T) 

+9x{^, Yo{t) + a7]{9rAto)){AT^ + /; /(e, yo(s) + <Tv{9s^Pe0o)) ds} , (13) 
Yi{0) = hl{0,u). 

Noticing that ri^{oj) and r]{-ip^uj) have identical distribution, together with the fact that 

1 Bs 1 /■* B(t-^.) 

r]'{uj) = ^ e~—dWs = ^ dWu = 7]%9tU}) , u = s + t, 

we see that r]'^[9tuj) and r]{i{ji;Uj) are identically distributed. Recall the transformation introduced 
in the beginning of this section X(t) = x{t) and Y(t) = y{t) — arj^{9tuj). We thus obtain a lower 
dimensional, slow stochastic system on the slow manifold, 

X = Ax + f{x,}f{x,9tuj) + ar]{ipeOj)), (14) 

Moreover, h'^{(^,oj) = h^{^,uj) +e/if(^,a;) is an approximation or a first order truncation of 
h^{^,uj), and hence we have an approximate slow system 

X = Ax + f{x,h'^{x,9tOj) + ar]{^plrUJ)). (15) 

This is the slow system we will work on for parameter estimation in the next section. 
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3 Parameter estimation on a random slow manifold 



If the original slow-fast system (l)-(2) contains an unknown system parameter a, but only the slow 
component x is observable, we can conduct parameter estimation using the slow system (15). In 
fact, this unknown parameter a is carried over to the slow system (15) which is now rewritten as 

X = Ax + f{x,k{x,etu}) + ar]{'4}e0j),a), x G M", (16) 

for £ sufficiently small. 

Since this slow system is lower dimensional than the original system, this method offers an ad- 
vantage in computational cost, in addition to a benefit of using only observations on slow variables. 
It is often more feasible to observe slow variables than fast variables [3]. 

Assume that we have observation, x^i,, on the slow component x only, and let us estimate the 
system parameter a. An estimator a^{e) for a may be obtained with existing techniques as reviewed 
in our earlier work [16] or [3, 6] and references therein. Here the subscript S indicates that the 
parameter estimation is conducted on the slow system (16). For example, a^(e) may be obtained 
by minimizing the objective function F{a) = E|| 

We then compare this estimator a|.(e) with the estimator, aE(e) (without superscript S*), based 
on the original slow-fast system (l)-(2), when observations Xq^, yoj, are available for both components 
x,y. For example, aE{e) may be obtained by minimizing the objective function I?'(Gf-) — ]E[||xo5 — 

A stochastic Nelder-Mead method is used to minimize the objective functions. In the next 
section, we demonstrate this method with an example. 



4 An example 

In this section, we demonstrate our parameter estimation method based on random slow manifolds 
by a simple example. 

Example 1. Consider a slow- fast stochastic system 

X = O.OOlx - axy, x(0) = xq G M, (17) 

y = -(-2^ + ?^^') + ^^*' y(0) = yoGM, (18) 
e dUU ^/e 

where a is a real unknown positive parameter, e is a small positive scale separation constant, a is 
a constant noise intensity, and Wt is a scalar Wiener process. 

See Figure 1 for a phase portrait of the corresponding deterministic system (u = 0). 

In this system, the nonlinear terms are not global Lipschitz. But if additionally it has an 
absorbing set, we can cut-off the nonlinearities without affecting the long time, slow dynamics 
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Figure 1: Deterministic dynamics - Phase portrait for x = O.OOlx — xy, y = ^{—y + ^q^"^) with 
e = |: The global attr actor is clearly seen (Red or thick curve within —7<x<7 and near the 
x— axis). 



(almost surely). Indeed, for arbitrary constants M and K, we have 

dx"^ = 2xdx = (0.002x2 _ 2ax^y) dt, 
d{My - Kf = 2{My - K)M dy + d[y, y] 

, 2M2 5 M2 „ 1MK MK 9 , , , a 

= + ^x^y + y - + )dt + 2M My -K) — 

^ e 300e e 300e e ' 

Therefore, 



d[^/^2a{My-Kf) 

1,M2 2 ^^^2, , 2aM2 2, 0.002M2 - 2aMif + MVe 2, 

= — ( x^ + 2a(My - Kf) dt y^ dt + —x^dt 

e^SOOe ^ e 300e 

_^2aK^ + 2aM^ + 4aM(M, - K)^ dW, (19) 

Taking M = e, K = then we see that 

0.002M2 - 2aMK + JVP/e 0.002^2 - 2 + e 

= < 0, tor small e, 

300e 300e 
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and 

= — ^ + 2ae(T 

£ ae-^ 

Thus 



|e(^x2 + 2aM\y - K/Mf) < -^^(^x^ + 2aM\y - K/Mf) + A + 2aea\ (20) 

By the Gronwah inequahty, we conclude that 

^^x^ + 2ae\y-^f) 
300 ae'^ ' 



^300 " ae^' ^ ae^ 

+ 2ae'(yo - \?y~-^ + (—2 
300 ae^ ae^ 



^ (w^o + 2ae2(yo-::;2 + (772 +2«^'^')- (21) 



This means, for fixed e and a > 0, the dynamics of the system (17)-(18) wih eventually stay in an 
ellipse (almost surely), i.e. there is a random absorbing set. We can then cut-off the nonlinearities 
outside this absorbing set to obtain a modified system which has, almost surely, the same long time, 
slow dynamics as the original system [7]. In the following calculations, we actually have omitted 
this cut-off procedure for simplicity. 

By the random transformation (4), SDEs system (17)-(18) are converted into the following 
system 

X = ^.m\X - aXiJ ^arf(et^)), (22) 

Y = --y + -— X^. (23) 
e £600 ^ ' 

Therefore, there exists an satisfying 

whose graph is a random slow manifold for the random system (22)-(23). In fact, h^{S^,u}) has an 
approximation /i^(^,tj) (with error C'(£^)), 

h^^uj) = — + e(-—Omi^ ^ a-—aa se'dWs). (25) 

' 600 ^ 300 180000 300 J _^ ' ^ ' 
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So the approximated slow system is 

X = Q.OOlx — ax (^ar]{tpi;Uj) + h'^{x,9tuj)). 



(26) 



with rj^tp^u}) = J_^e^ dWs{tpe^)- We will illustrate that the parameter estimator for a based on 
this low dimensional, reduced system (26) is a good approximation of the parameter estimator 
based on the original system (17)-(18). 

Slow-fast system and random slow manifold The random slow manifold is the graph of 
h^{^,oj) = arj{Tp!.oj) + h^{^,uj), where is as in (24). It is a curve (depending on random samples). 
The random orbits of the system (17)-(18) approaches to this curve exponentially fast. Figures 
2-3 show some orbits of the slow-fast system (17)-(18) and the approximate random slow manifold 
f)^(^,a;) = ar]{tp!rU}) + /i^(^,cj), where is as in (25), with different a and e values. 

The orbits of the system (17)-(18) decay quickly to the random slow manifold. 

Figure 4 shows several samples or realizations of the random slow manifold with different a 
values. 

Nelder-Mead method and stochastic Nelder-Mead method The deterministic Nelder- 
Mead method (NM) is a geometric search method to find a minimizer of an objective function 
F{a). Starting from an initial guess point, it generates a new point (reflection point, expand point, 
inside/outside contraction point or shrink point) by comparing function values, and thus get a 
better and better estimator until the smallest objective function value in this iteration reaches 
the termination value (prescribed error tolerance). The algorithm of the NM method and its 
improvements have been widely studied and utilized [1, 9, 10]. The method has its advantage that 
the objective function need not to be differentiable, and it can thus be used in various applications. 
As noted in [ 1, I I], the Nelder-Mead method is a widely used heuristic algorithm. Only very limited 
convergence results exist for a class of low dimensional (one or two dimensions) problems, such as 
in the case when the objective function F is strictly convex [S]. For the numerical simulation, with 
analytical objective functions, Matlab function fminsearch can be used to find a minimizer. 

However, when dealing with problems with noise, NM method has the disadvantage [2, 1] that 
it lacks an effective sample size scheme for controlling noise, as shrinking steps are sensitive to the 
noise in the objective function values and then may lead to the search in a wrong direction. In fact 
an analytical and empirical evidence is known for the false convergence [2] on stochastic function. 
So we use the stochastic Nelder-Mead simplex method (SNM) [4] to mitigate the possible mistakes 
in the stochastic setting. The newly developed Adaptive Random Search in [ I] consists of a local 

. Nik) 

search and a global search. It generates a new point and new objective function F = ^ Fi/N{k) 

i=l 

in the A;— th iteration with increasing number of the sample size scheme N(k). A proper choice 

N{k) 

for N{k) is [^/k], with [c] the largest integer not bigger than c. Here ^ Fi is the sum of N{k) 

i=l 
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objective function values for F. SNM leads to the convergence of F at search points to minF(a) 

a 

(and thus we obtain a minimizer a*), with probability one. 

Parameter estimation To illustrate our method for parameter estimation on the random slow 
manifold, we fix e = 0.01 and a = 0.01 in the following numerical experiments. 

We want to estimate the parameter a by using both the original slow-fast system and the slow 
system, in order to demonstrate that the slow system is appropriate for parameter estimation, when 
e is sufficiently small. 

Step 1: Generate observations 

Take the true value a = 0.1 (say) and numerically solve (17)-(18) with an initial condition 
{xq, uo) to get J samples of observational data (x^"'^, j = 1, • • • , J at time instants ti, i = 1, ■ ■ ■ , I 
(save these data). 

Step 2: Estimator for the original slow- fast system 

Take two initial guesses for the unknown system parameter a = and a = oi randomly and 
solve the original slow-fast system (17)-(18), with the same initial condition (xq, uq) and time points 
ti, i = 1, ■ ■ ■ , F Thus we obtain and values which depend on a. 

Using the stochastic Nelder-Mead Algorithm, we find a parameter estimator qe such that the 

I J . .. 

objective function F{a) = E ^ ^ ((x* — x^],)^ -|- (y* — y*],)^) is minimized. 

i=ij=i 

Step 3: Estimator for the slow system 

Take two initial guesses oq and ai randomly and solve the slow system (26) with the same initial 
condition xq at the same time instants U, i = 1, ■ ■ ■ ,F We thus obtain X* which depends on ao or 
ai. 

By the stochastic Nelder-Mead method as in Step 2, we find the parameter estimator by 

I J , 

minimizing the objective function F(a) = E ^ ^ (X* — x^-'^)^. 

i=ij=i 

Numerical experiments Figure 5 shows the objective functions F{a) and F(a) with true a = 0.1 
(left) and a = 1 (right). Here we used 30 paths to get the expectation in the objective function. 

For £ = 0.01 and a = 0.01, Figures 6 and 7 are objective function values and estimators of each 
iteration for slow-fast system (17)-(18) (top) and reduced system (26) (bottom) with a = 0.1 and 
a = 1, respectively. 

We observe that, as proved in [4], the objective function value F{a) tends to minF(a)(=0), 

a 

while the minimizer or the parameter estimator qe provides an accurate estimation of the system 
parameter a. For sufficiently small e > 0, the objective function value F(a) for the slow system 
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will also get closer and closer to 0, and the minimizer or the parameter estimator is a good 
approximation of ge- 
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Figure 2: Orbits of slow-fast system (17)-(18) (blue curves) and its slow manifold expansion 
f)^(^,a;) = ari{ip!rUj) + h^{S^,uj) (red curve), where is as in (25), with a = 0.01 and e = 0.01: 
a = 0.1 (left ) and a = 1 (right). 
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Figure 3: Orbits of slow-fast system (17)-(18) (blue curves) and its slow manifold expansion 
f)'^(^,a;) = ar]{ip!rUj) + h^{^,uj) (red curve), where fi^ is as in (25), with a = 0.05 and e = 0.01: 
a = 0.1 (left ) and a = 1 (right). 
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Figure 4: Random slow manifold expansion f)''(^,a;) = (Tr]{il)i,uj) + h^{(^,uj), where is as in (25) 
with e = 0.01: a = 0.01 ( left ) and a = 0.05 (right). 
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Iteration Iteration 

Figure 6: Objective function F{a) and estimator aE (top) for slow-fast system (17)-(18); and 
objective function F(a) and estimator (bottom) for slow manifold reduced system (26): a = 0.01, 
e = 0.01 and true value a = 0.1. 
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Iteration Iteration 

Figure 7: Objective function F{a) and estimator ue (top) for slow-fast system (17)-(18), and 
objective function F(a) and estimator (bottom) for slow manifold reduced system (26): a = 0.01, 
e = 0.01 and true value a = 1. 
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